CMRSegTools: An open-source software enabling reproducible research in segmentation of acute myocardial infarct in CMR images

In the last decade, a large number of clinical trials have been deployed using Cardiac Magnetic Resonance (CMR) to evaluate cardioprotective strategies aiming at reducing the irreversible myocardial damage at the time of reperfusion. In these studies, segmentation and quantification of myocardial infarct lesion are often performed with a commercial software or an in-house closed-source code development thus creating a barrier for reproducible research. This paper introduces CMRSegTools: an open-source application software designed for the segmentation and quantification of myocardial infarct lesion enabling full access to state-of-the-art segmentation methods and parameters, easy integration of new algorithms and standardised results sharing. This post-processing tool has been implemented as a plug-in for the OsiriX/Horos DICOM viewer leveraging its database management functionalities and user interaction features to provide a bespoke tool for the analysis of cardiac MR images on large clinical cohorts. CMRSegTools includes, among others, user-assisted segmentation of the left-ventricle, semi- and automatic lesion segmentation methods, advanced statistical analysis and visualisation based on the American Heart Association 17-segment model. New segmentation methods can be integrated into the plug-in by developing components based on image processing and visualisation libraries such as ITK and VTK in C++ programming language. CMRSegTools allows the creation of training and testing data sets (labeled features such as lesion, microvascular obstruction and remote ROI) for supervised Machine Learning methods, and enables the comparative assessment of lesion segmentation methods via a single and integrated platform. The plug-in has been successfully used by several CMR imaging studies.


Introduction
This appendix provides a brief description of the computational methods implemented within CMRSegTools for the segmentation of lesions in LGE MR images.

Methods
All the methods described below are performed within a subset of the image pixels selected by a prior LV myocardium segmentation method. The LV segmentation returns a non-empty set of pixels inside the endocardial and epicardial contours, where p denotes a pixel at (i, j) image coordinates.
Working on a subset of the image pixels reduces the computational complexity O(N ) (expressed as a function of the image size N ) as the algorithm only computes the solutions on the pixels of interest.

Manual cut-off (n-SD)
This is the most used semi-automatic method based on STRM. Several studies have demonstrated good correlation with pathology. This method requires the definition of a ROI in the remote healthy myocardium (figure 1a). Therefore, the threshold is calculated by where µ is the mean of the signal intensity in the remote healthy myocardium and σ is the standard deviation, c is a positive integer usually between 2 and 10. The subset of the lesion pixels are selected as where I(p) is the signal intensity of pixel p. Endocardium (red) contours, the LV/RV junction landmark (yellow), the lesion segmentation (magenta) and the remote healthy myocardium (turquoise); (b) The CMRSegTools GUI displays the histogram highlighting the bins from the remote ROI (turquoise) and lesion segmentation (magenta).

Manual histogram-based segmentation
All the pixels which have an intensity value within the range {a, b} are categorised as lesion where a and b are the range limits set by the user on the GUI with the lower and upper cursors on the histogram or using the Histogram Range text boxes (2b).

Full-Width at Half-Maximum (FWHM)
In this automatic method, the threshold is calculated as where I max is the maximum signal intensity within the myocardium (LV seg ). CMRSegtools includes the following options: • FWHM Max. The threshold is calculated from the maximum SI within the myocardium of the k-th slice in K slices.
• FWHM Region. This method requires the definition of a ROI in the healthy myocardium. The threshold is calculated as where µ is the mean of the signal intensity in the remote myocardium.
• FWHM Region3D. The method is a generalisation of the FWHM Region. The threshold is calculated with the maximum SI in the k-th slice and the mean SI of a unique remote ROI defined in slice k remote .

Gaussian Mixture Model
This approach models the histogram of intensity values x (i.e. SI) in the myocardium as where w R corresponds to the relative proportion of healthy tissue and w G to the lesion respectively (Rician-Gaussian mixture model [1]). The probability density function (PDF) of the Gauss distribution is where µ and σ 2 are the mean and variance of the distribution respectively, parameters to be identified by the algorithm. The PDF of the Rice distribution is where m is a parameter and I 0 (z) is the modified Bessel function of the first kind with order zero [2]. The algorithm relies on an Expectation-Maximisation (EM) method to estimate the parameters of the mixture model.

Feature Analysis and Combined Thresholding
The feature analysis and combined thresholding (FACT) algorithm proposed by Hsu et al. in [3] is a workflow composed of the following steps: 1. Histogram clustering. The GMM method (Section 2.4) is applied in order to find the mean µ and the standard deviation σ of the healthy myocardium and define a binary image by applying where I(k, p) is the signal intensity of pixel p at slice k, and T = µ + cσ.

2.
Region-Based Feature Analysis. On each slice k, isolated 2D regions are grouped to 3D volumes (stack of 2D slices) by using a neighbour connectivity criteria providing a 3D binary image (volume) with the voxels of interest. In order to remove false positive voxels, the following filters are performed: • Volume filter removes regions smaller than a minimum setting based on the myocardial mass density.
• Distance filter removes any region more than 2 mm away from the endocardial border.
• Intensity filter removes voxels about 50% darker than the average intensity of a 3D region.
3. FWHM Thresholding. In order to remove voxels from partial volume tissues, a threshold defined as T = Imax+µremote 2 is applied. Steps 2 and 3 can be repeated in order to reduce false positives voxels/regions.

4.
Detection of the MVO regions. A morphology closing filter is used to identify small holes surrounded by bright infarcted voxels. These small holes are then labeled as no-reflow regions.

Hidden Markov Random Field model with Expectation-Maximisation (HMRF-EM)
This method has been developed based on the algorithm proposed by Zhang et al. [4] in which the objective is to define the corresponding labels x = {x 1 , x 2 , x 3 , ..., x N } for the pixels in an image y = {y 1 , y 2 , y 3 , ..., y N } where y i is the SI of the i-th pixel. In the case of CMR, the set of labels L in the myocardium (LV seg ) is defined as L = {l 1 , l 2 , l 3 , ..., l M } (namely L = {normal, lesion}). Within a Maximum likelihood framework, the method relies on the solution to the following optimisation problem: where Θ is the set of parameters (mean and standard deviation) of the PDF that defines each label; which is written formally as Θ = {θ l = (µ l , σ l ) | l ∈ L}. The estimation of these set of parameters is made using the EM algorithm.
The HMRF-EM workflow in CMRSegTools includes the following steps [5]: 1. Cartesian to polar transformation of the LV seg region. A polar space is better adapted to the myocardium morphology. In the polar space neighbour operations are more effective as pixels in the same circumference (cartesian space) are aligned in the same line (polar space). Figure 4 shows the output of the Cartesian to polar transformation.
2. Initial classification. An initial label map and the corresponding parameters (mean and standard deviation) are computed by finding the SI which is most likely to be lesion. This initial classification uses a threshold defined as T = 4×Imax 5 with I max given by 4. This step provides the initial parameter set Θ (0) for the label map x (0) . 3. EM algorithm. The t-th iteration of the process, the Expectation step determines the conditional expectation of the parameters and the Maximisation step calculates the next estimation of parameters Maximum a posteriori (MAP) estimation is performed for calculating the labels based on the current set of parameters Θ (t) where U (y | x, Θ (t) ) is the likelihood energy distribution. The EM process continues until the energy is stable or a maximum number of iterations is reached ( Figure 5). 4. Post-processing and identification of MVO regions. The LV seg region is transformed from polar to Cartesian. A connected component analysis is applied to remove small clusters (i.e. less than 0.1 g) or isolated clusters too distant from the endocardium (distance above 1.5 mm), removing false positives. As the MVO region can be completely inside the lesion region or at the edge of myocardial boundaries, connected components which size is over 50% of the myocardium size are removed. External bordering pixels that are not labeled as lesion are removed if they are above a defined threshold (set to 50%). Figure 6 illustrates the HMRF-EM workflow. Figure 6: HMRF-EM workflow [5].